home *** CD-ROM | disk | FTP | other *** search
/ Delphi Programmer's Power Pack / Delphi Volume 1.iso / s_to_z / statone / stat.pas < prev    next >
Pascal/Delphi Source File  |  1996-09-15  |  5KB  |  201 lines
  1. unit Stat;
  2.  
  3. interface
  4.  
  5. uses
  6.   SysUtils,
  7.   WinTypes,
  8.   WinProcs,
  9.   Messages,
  10.   Classes,
  11.   Graphics,
  12.   Controls,
  13.   Forms,
  14.   Dialogs,
  15.   dynary;
  16.  
  17. const
  18.     STATVERSION = 'Ver 1.10';
  19.   MISSING = 5e-324;  {not used in this version}
  20.  
  21.   type
  22.   TStatResults = record
  23.     N, NMissing: LongInt;
  24.     Mean, Median, Sum, SS, Variance, Skewness, Kurtosis, StdErr,
  25.     StdDev, GeoMean, HarMean, Bisquare, IQR: Double;
  26.  
  27.   end;
  28.  
  29.   StatOne = class(TDoubleArray)
  30.   private
  31.     { Private declarations }
  32.     FAboutStatOne: string;
  33.     FStatResults: TStatResults;
  34.     procedure SetAboutStatOne(value: string);
  35.  
  36.   protected
  37.     { Protected declarations }
  38.   public
  39.     { Public declarations }
  40.       property StatResults: TStatResults read FStatResults;
  41.     procedure CalcStat;
  42.     procedure Bisquare(const c: integer);
  43.     function Percentile(Percent: double): double;
  44.  
  45.  published
  46.     { Published declarations }
  47.     constructor Create(AOwner: TComponent); override;
  48.     property AboutStat: string read FAboutStatOne write SetAboutStatOne;
  49.  end;
  50.  
  51. procedure Register;
  52.  
  53. implementation
  54.  
  55. procedure StatOne.Bisquare(const c: integer);
  56.     function MaxAbs(n1, n2: double): double;
  57.   begin
  58.         if n1 = 0.0 then begin
  59.       Result := Abs(n2);
  60.       Exit
  61.       end else if n2 = 0.0 then begin
  62.       Result := Abs(n1);
  63.       Exit
  64.       end;
  65.     if Abs(n1) >= Abs(n2) then
  66.         Result := Abs(n1)
  67.     else
  68.         Result := Abs(n2);
  69.   end;
  70.  
  71. var
  72.     b1, b2, s, u, w, wx, weight: double;
  73.   i: integer;
  74. begin
  75.     CalcStat;
  76.     b2 := FStatResults.Mean;
  77.     s := c * FStatResults.IQR / 2;
  78.   while (Abs(b1 - b2) / MaxAbs(b1, b2)) > 0.01 do begin
  79.     b1 := b2; w := 0; wx := 0;
  80.     for i := 1 to FStatResults.N do begin
  81.       u := (Value[i] - b1) / s;
  82.       if abs(u) <= 1 then begin
  83.         weight := sqr((1 - sqr(u)));
  84.         w := w + weight;
  85.         wx := wx + (weight * Value[i])
  86.       end;
  87.     end;
  88.     b2 := wx / w;
  89.   end;
  90.   FStatResults.Bisquare := b2;
  91. end;
  92.  
  93. procedure StatOne.CalcStat;
  94. var
  95.     i: integer;
  96.   N, NGeo: LongInt;
  97.   x, ssLn, ssInv, v, m1, m2, m3, m4, Temp: Double;
  98.  
  99. begin
  100.     FStatResults.N := 0;
  101.     FStatResults.NMissing := 0;
  102.   FStatResults.Sum := 0;
  103.   FStatResults.SS := 0.0;
  104.   FStatResults.Skewness := 0.0;
  105.   FStatResults.Kurtosis := 0.0;
  106.   FStatResults.StdErr := 0.0;
  107.   ssLn := 0.0; ssInv := 0.0;
  108.   v := 0.0; m1 := 0.0; m2 := 0.0; m3 := 0.0; m4 := 0.0; N := 0; NGeo := 0;
  109.  
  110.   Sort;
  111.     for i := 1 to Size do begin
  112.       if i > MISSING then begin
  113.         Inc(N);
  114.       v := (Value[i] - m1) / N;
  115.       m4 := m4 - (4 * v * m3) + (6 * v * v * m2) + ((N * N) - (3 * (N - 1))) * N * (N -1) * v * v * v * v;
  116.       m3 := m3 - (3 * v * m2) + (N * (N - 1) * (N - 2) * v * v * v);
  117.       m2 := m2 + (N * (N - 1) * v * v);
  118.       m1 := m1 + v;
  119.       FStatResults.Sum := FStatResults.Sum + Value[i];
  120.       FStatResults.SS := FStatResults.SS + Sqr(Value[i]);
  121.       if Value[i] > 0 then
  122.       begin
  123.           ssInv := ssInv + 1 / Value[i];
  124.           ssLn := ssLn + Ln(Value[i]);
  125.         Inc(NGeo);
  126.       end;
  127.     end else
  128.          Inc(FStatResults.NMissing);
  129.   end;
  130.       FStatResults.N := N;
  131.     FStatResults.Mean := m1;
  132.     if N > 1 then
  133.             FStatResults.Variance := m2 / (N - 1)
  134.     else
  135.       FStatResults.Variance := 0.0;
  136.     FStatResults.StdDev := sqrt(FStatResults.Variance);
  137.     FStatResults.Sum := m1 * N;
  138.  
  139.     if N > 2 then
  140.         FStatResults.Skewness := (N * m3) / ((N - 1) * (N - 2) *(FStatResults.StdDev * FStatResults.StdDev * FStatResults.StdDev))
  141.     else
  142.       FStatResults.Skewness := 0.0;
  143.  
  144.     if N > 3 then
  145.     begin
  146.         Temp := (N - 1) * (N - 2);
  147.         FStatResults.Kurtosis := ((N * (N + 1)* m4) - (3 * m2 * m2 * (N - 1))) /
  148.           (Temp * (N - 3) * FStatResults.StdDev * FStatResults.StdDev *
  149.         FStatResults.StdDev * FStatResults.StdDev)
  150.     end
  151.     else
  152.     begin
  153.       FStatResults.Kurtosis := 0
  154.     end;
  155.  
  156.     FStatResults.StdErr := FStatResults.StdDev / sqrt(N);
  157.     FStatResults.GeoMean := Exp(ssLn / NGeo);
  158.     FStatResults.HarMean := NGeo / ssInv;
  159.     FStatResults.Median := Percentile(50);
  160.     FStatResults.IQR := Percentile(75) - Percentile(25);
  161.   end;
  162.  
  163.   function StatOne.Percentile(Percent: double): double;
  164.   var
  165.       x: double;
  166.     i, N: LongInt;
  167.   begin
  168.       if (Percent < 0.0) or (Percent > 100.0) then begin
  169.         showmessage('Percent out of range');
  170.         exit;
  171.     end;
  172.     N := FStatResults.N;
  173.       Percent := Percent / 100.0;
  174.     x := Int(N * Percent);
  175.     if x = (N  * Percent) then begin
  176.       i := Trunc(x);
  177.       Result := (Value[i] + Value[i+1]) / 2;
  178.     end else begin
  179.       i := Round(N * Percent) {+ 1};
  180.       Result := Value[i];
  181.     end;
  182.   end;
  183.  
  184. constructor StatOne.Create(AOwner: TComponent);
  185. begin
  186.   inherited Create(AOwner);
  187.   FAboutStatOne := STATVERSION;
  188. end;
  189. procedure StatOne.SetAboutStatOne(value: string);
  190. begin
  191.     FAboutStatOne := STATVERSION;
  192. end;
  193.  
  194.  
  195. procedure Register;
  196. begin
  197.   RegisterComponents('Ted', [StatOne]);
  198. end;
  199.  
  200. end.
  201.